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■ Abstract 

Oh, We study the phase structure of the random-plaquette Z2 lattice gauge model in three dimensions. In 

I . this model, the "gauge coupling" for each plaquette is a quenched random variable that takes the value 

^ ■ /3 with the probability 1 — p and — /3 with the probability p. This model is relevant for the recently 



proposed quantum memory of toric code. The parameter p is the concentration of the plaquettes 
with "wrong-sign" couplings — /3, and interpreted as the error probability per qubit in quantum code. 
In the gauge system with p = 0, i.e., with the uniform gauge couplings /3, it is known that there exists 
a second-order phase transition at a certain critical "temperature", r(= /3^^) = Tc = 1.31, which 
separates an ordered (Higgs) phase at T < Tc and a disordered(confinement) phase at T > T^- As -p 
increases, the critical temperature T^ij)) decreases. In the p — T plane, the curve Tc(p) intersects with 
the Nishimori line Tn(p) at the certain point (pc, Tn(pc)). The value pc is just the accuracy threshold 
for a fault-tolerant quantum memory and associated quantum computations. By the Monte-Carlo 
simulations, we calculate the specific heat and the expectation values of the Wilson loop to obtain 
the phase-transition line T(.(j>) numerically. The accuracy threshold is estimated as Pc — 0.033. 
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1 Introduction 



Gauge theory plays a very important role in modern physics. All the interactions between elementary 
particles are described by gauge theories. Furthermore, gauge theories present important and clear 
understandings for other areas of physics including superconductivity, fractional quantum Hall effect, 
random spin systems, localization by random vector potentials. 

Recent dcvelopements of nanotechnology elevate the old interesting theoretical ideas to realistic 
objects. Quantum computer(QC) is one of them. One of the most difficult problem for constructing 
QC is how to keep quantum states stable against noises while quantum computations are performed, 
i.e., the problem of decoherence. In order to make a quantum computation fault-tolerant, the 
involved qauntum states should have certain "quantum numbers" , and these quantum numbers are 
required to stay as good ones during the computation. From this point of view, one may expect 
that gauge-theory mechanisms like Aharonov-Bohm(AB) effect play an important role in quantum 
computation. 

In fact, KitaevpP recently proposed a fault-tolerant quantum memory and quantum computations 
that are based on the AB effect of discrete Z2 gauge symmtery. After that many interesting studies 
on Kitaev's model appeared[21ElSElEllZllHl- Among them, two of the present authors 7| generalized 
Kitaev's model and also clarified the relationship between Kitaev's model and a spontaneously broken 
U{1) gauge theory. 

For quantum computations, it is important to estimate its accuracy threshold. As discussed in 
Ref.0], it is obtained for Kitaev's model by studying the random-plaquette gauge model(RPGM) 
on a three-dimensional {3D) lattice with Z2 gauge symmetry. In the RPGM, the (inverse) gauge 
coupling for each plaquette takes the values ±/3 with random sign; the probability to take (3 is given 
by 1 — p and the probability of "wrong-sign" — /3 is p. The original nonrandom {p — 0) 3-D Zi 
lattice gauge theory is known to be dual to the classical 313 Ising model. The fluctuations of gauge 
field are controlled by the "temperature" (the coupling constant) T = f3^^. The system exhibits a 
second-order phase transition at the critical temperature T = Tc — 1.31. For T < Tc the gauge- 
field fluctuations are small and the system is in the ordered Higgs phase, whille for Tc < T, the 
fluctuations are large and the system is in the disordered confinement phase _9j. As p increases from 
p ~ 0, the critical temperature Tdp) decreases due to the randomness in gauge couplings, and it 
shall vanish at a certain value pa, Tc{pq) — 0. In Fig.l, we illustrate a schematic phase diagram 
in the p — T plane. The solid line expresses Tc{p), which may decrease very rapidly near pq. The 
dashed line is the Nishimori line T{p) — Tn(p)|10i on which fluctuations of gauge field caused by the 
quantum effect and those caused by the random couplings are equal. As we explain it in detail later 
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on, we consider a particular method of error corrections in QC's that is reahzed along the Nishiniori 
line. Also the system should be in the ordered Higgs phase in order that the quantum computations 
are fault-torelant. Thus the intersection point of the critical line Tdp) and the Nishimori line, i.e., 
Tc{Pc) = 7n(pc)j locates the accuracy threshold for the quantum computations by Kitaev's model. 
One concludes that the quantum computations are fault-torelant for p < pc- 

The present paper is organized as follows. In Sec. 2, we shall review Kitaev's model of a quantum 
toric code and its relation to the Z2 gauge model. In Sec. 3, the relationship between the accuracy 
threshold of the toric code and the random spin and gauge models is discussed. Section 4 is the 
main part of the present paper. We present the results of Monte-Calro simulations of the RPGM. 
We calculate the specific heat and expectation values of the Wilson loop for various points in the 
p ~ T plane. There are (at least) two phases in the model; one is the ordered Higgs phase and the 
other is the disordered confinement phase. From these studies, we estimate the accuracy threshold 
of the toric code as Pc — 0.033. This value of threshold is consistent with the previous studies^ 
In Ref. 4 , the Z2 RPGM in 3-D was investigated by studying the lowest-energy configurations for 
various random gauge couplings, and then the critical concentration along the T — Q axis was 
estimated as Pq = 0.029. By using specific symmetry of the model on the Nishimori line, it is proved 
that the above value gives a lower bound for the accuracy threshold pc 10 . On the other hand in 
Ref.0, the four- dimensional RPGM was studied. The value of p at the intersection point of the 
phase boundary and the Nishimori line of this model gives an upper bound of pc of the RPGM in 
3-D. By using the self duality of the model, it is obtained as 0.11, which is much larger than our 
result of Pc- Section 4 is devoted for conclusion. 

2 Quantum toric code 

The quantum toric code was first proposed by Kitaev^, and many interesting studies on that 
idea appeared after that. It is more transparent to formulate Kitaev's model as a Z2 lattice gauge 
theory 7 . The quantum toric code is viewed as a set of qubits put on links of a two-dimensional(2-D) 
square lattice with the periodic boundary condition (i.e., a torus). Each link of this lattice is denoted 
as (x, i) where x is the site index and i = 1, 2 is the spatial direction index. The qubit residing on the 
link [x, i) is described by a quantum Z2 variable Z^i and its conjugate variable X^i- Z^i is regarded 
as the Z2 gauge field, and so X^i plays the role of the electric field. They satisfy the commutation 
relations like 

XxiZxi = —ZxiXxi, XxiZyj = ZyjXxi for {x,i) ^ {y,j)- (1) 



3 



In the matrix representation, Z^i = <7%^ and X^i — cri^j, where ct^ and cr^ are 2x2 Pauh spin 
matrices. 

Let us consider the following generalized Hamiltonian, 
Br = Hz + Ht + H'^z^ 

Hz = -Ai^X„-A2 Yl ZZZZ + Y^.c, 

link plaqucttc 

= -lY,^^l+.,Z,^i^. + MY,i;li^, + \^.c., 

link site 

H^z = -7^<^l+.^..<^.+M^<^I<^.+H.c., (2) 

link site 

where Ai and A2 are coupling constants, and ijjx and (fx are "relativistic" fermions sitting on site a; 
with the mass M and the hopping parameter 7. They satify the canonical anticommutation relations 
like 

{■4'x,i^l] = Sxy, Wx,fl} = Sxy (3) 

The above Hamiltonian ^ is invariant under Z2 local gauge transformation. In fact, we have 

{y,i)ex 

[Gx,Ht]^0, (4) 

where (y, i) d x denotes 4 links emanating from site x, and the operator Gx is the generator of gauge 
transformations. Respecting the gauge invariance of Hx, the phyical states \phys) should be gauge 
invariant, i.e., satisfy the following physical-state condition, 

Gx\phys) = \phys). (5) 

The Hamiltonian of Kitaev's model is obtained by setting Ai = 0, A2 = 1,7 = and M = 2[7], 
and then the quantum space of the gauge field and that of the fermions are coupled only through 
the gauge-invariance condition. In this case, it is easy to see that the ground states \GS)k satisfy 
the following conditions, 

V'iV.|G5)K-0, ^lipx\GS)K = 0, 

Ax\GS)k = \GS)k, Bp\GS)k^\GS)k, (6) 
for all the sites x and plaquettes P, where the operators and Bp are defined as follows, 

= n ^y^' ^^'^ = n ^v^- 

{y,i)ex (y,i)eP 
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These operators are called stabilizers or check operators. From Eq. Q , is the generator of a local 
gauge transformation of the Z2 gauge field Z^i, and Bp measures Z2-gauge flux penetrating the 
plaquette P, i.e., Bp = 1 for fluxless states and Bp = —1 for fluxful states. Eq.® means that 
\GS)k involve no fermions (pure-gauge state), gauge-invariant, and fluxless. 

As the system is defined on the torus, the ground states are not unique but 4-fold degenerate. 
They satisfy Eq.® and are distinguished by a pair of "nonlocal" quantum numbers Za{a — 1,2) 
defined as 

Za=Y[Z,,, (8) 

where Cf(a — 1,2) are two noncontractible closed loops on the original lattice (see Fig. 2). The 
ground states have the quantum numbers {Zi, Z2) — (1, 1), (1, —1), (~1, 1) and (—1, —1). Quantum 
informations for the toric code are stored in the above 4-fold degenerate ground states which work as a 
fault-tolerant qudit. Instead of Za, we can consider a set of conjugate quantum numbers Xa{a = 1,2) 
given by 

Xa=llX,,, (9) 

where Cx are nontrivial closed loops on the dual lattice (see Fig. 2) and X^i's in Q reside on the links 
of the original lattice that are crossed by Cx- It is easily verified that Za and Xa{a = 1,2) satisfy 
the same commutation relations with the Pauli matrices. Then the ground states are characterized 
by the quantum numbers {Xi, X2) — (1, 1), (1, —1), (—1, 1) and (—1, —1). 

The gauge-theory aspect of Kitaev's model becomes clear when we consider excitations in the 
system. Excitations obviously break the local condition ® at some specific sites and/or plaquettes. 
As the stabilizers satisfy the identities like 

Y[ A, = 1, and II Bp = l, (10) 

all sites all pi's 

these excitations should appear in pairs. A simple example is the fermion-pair state with two 
fermions at x and y, 

\F;C,y)^^Pl[l[z)^l\GS)K, (11) 

where C^y is a certain path on the original lattice connecting x and y. Another example is the 
vortex-pair state with two vortices at dual sites x* and y* , 

\V;C,'y^)^{ n x)\GS)k, (12) 

where Cx*y» is a certain path on the dual lattice connecting x* and y*, and X's in 1)12(1 are on 
the links crossing Cx'y-{see Fig. 3). Wave function of the state containing the above fermionic and 
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vortex excitations acquires nontrivial phase factors when one of the excitations winds around the 
other. This is the AB effect of the gauge theory. 

As explained above, quantum informations are stored in a set of the 4-fold degenerate ground- 
states, L = {|G5)}. During a quantum computation, quantum states in quantum memory are to 
be damaged by noises which cause decoherence of the quantum states. As a result, excited states 
appear, i.e., errors occur. One can check whether errors have occured by measuring the operators 
tpl;tpx, ftfx, and Bp in the quantum states and judging whether the conditions (jSJ are satisfied 
or not. We can correct these errors by applying suitable operators on the quantum states. For 
example, the state with vortex- type excitations can be corrected by applying the operator (n^)i 
a product of Xxi along the path appeared in the excited state, owing to the identity X'^^ = 1. On 
the other hand, for the charge-type excitations, the operator V'j^dl Z)lPx does the job. But we note 
that we can know only the locations of vortices and/or charges but not that of the error chain itself. 
Therefore we should consider specific methods of the error corrections. As a results there exists an 
accuracy threshold Pc, i.e., if the error probability per qubit and per unit computer time exceeds Pc, 
we cannot perform successful error corrections in practical manner. The threshold can be estimated 
by studying the random-bond Ising model(RBIM) in 2D and/or the RPGM in 3Z3, as we discuss it 
in the following section. 

3 Accuracy threshold and RPGM 

It is interesting and important to estimate the accuracy threshold Pc of the toric quantum code. The 
value of Pc should be compared with those of other quantum memories and quantum gates which 
are estimated as 10~^ to 1/300^]. In the most of the previous studies, concatenated quantum error 
correcting codes were considered instead of the surface codes in the present study. 

As discussed in detail in Ref . jSl El , Pc of the toric code can be obtained by studying the RBIM 
in 2D and/or the Z2 RPGM in iD. We note that pc depends on the practical methods of error 
corrections which one employs. In the method of Ref.^ it is estimated as pc = 0.029. Here we shall 
review the estimation of pc by using the statistical-mechanical models in two steps. In the first step, 
we assume that every measurement of check operators done in order to check the quantum states 
is perfect, i.e., without any errors. Then the relevant model to calculate Pc is shown to be the 2D 
RBIM. In the second step, we take into account the possible errors in measuring check operators 
themselves. In this case, the relevant model becomes the iD Z2 RPGM. The latter model resembles 
the 2D RPGM of the toric code discussed in the previous section, but the dimensions of two models 
are different. The third dimension is introduced as the time direction along which we measure check 
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operators successively. 

Let us start the first step. There are two types of errors; (i) phase errors associated with fermion- 
pair states of and (ii) spin-flip errors associated with vortex-pair states of (|12|l . Since these two 
types of errors can be discussed in the same manner, we focus here on the phase errors. A typical 
examlpe of phase-error syndrome is given in Fig. 4. The sites marked by filled circles denote error 
sites on which the check operator has the eigenvalue (—1), while = 1 on the remaining sites. 
This implies that a chain of the phase errors have occurred. In Fig. 5, we show one of the simplest case 
of two phase-error sites, where the error chain containing all error links is denoted by E. Therefore 
on the boundary sites of E, A^ = — 1. In order to correct these errors, we have to connect the error 
sites by a certain chain which we call E' , and apply a phase operator (O^^ri) along E' upon the 
quantum state under problem. There is of course ambiguity in choosing a suitable E' . Let us see 
this in detail. Because E plus E' form a closed loop, i.e., a cycle, the corrected state and the original 
state differ by the factor Hb+b' ^^i- this closed loop belongs to the homologically trivial class of 
the set of closed loops on the torus, we have 

n n Bp^i, (13) 

E+E' P with in E+E' 

beause of — 1 and the last equation of Eq. . Thus the correction is successful. On the other 
hand, if the closed loop belongs to the homologically nontrivial class, the first equality of H13|l does 
not hold, hence the above correction procedure does not work. 

Let us consider a toric code of an arbitrary large size and assume that the probability of error 
per qubit is p. An error chain E is characterized by a function of each link nE{t} {£ — link) like 
"-bI^) = 1 for ^ G £^ and nE{£) = for £ ^ E. As we explained above E' = E + C where C is a cycle 
(i.e., a chain with no boundary) and then we define a function nc{t) for C in the exactly same way 
with nE{(-)- The probability that error chain E occurs is given as 

Woh{E) = ]J(1 -p)i-"^(V*^^- (14) 

£ 

Let us consider the probability distribution for an arbitrary chain E' that has the same boundary 
with E, prob(£"|£')^. For £ e E' and i ^ E, the hnk functions are nci£) = 1 and n_E(^) =0 and its 
probability is given as 

(15) 



1-p 

On the other hand for £ ^ E' and £ ^ E with the link functions nc{£) = 1 and nE{£) = 1 

1-P 

P 



(16) 



^The actual probability for the chain E' is given by prob(i?) ■ prob(i?' 1_E). 
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For nc{i) = 0, the probability is unity. Then the conditional probability is given by (up to an 
irrelevant constant factor) 

proh{E) ■ pioh{E'\E) ^Y[cxp{j7jeui), (17) 

e 

where 

ue^l-2nc{i) e (18) 

f I, (.iE 

m = { (19) 
[ -1, i<^E, 

and e^^'^ —p/{l —p). 

From Ea. H18|l and the fact that C is a close chain, 

n^^=i' (20) 

at each site x. The above constraint 12U|) can be easily solved by introducing classical Ising spin 
variables a^' = ±1 at sites of the dual lattice, 

M£ — CTx-'CTx-'+i, (21) 

where the notation is obvious. Then the generating function of the probability distribution of the 
error chains E' is given by the following partition function Z of the RBIM, 

T]£ax'Crx*+i)- (22) 

The relation e^'^'^ — p/(l ^ p) between the temperature and the concentration of the "wrong-sign" 
bonds p defines the Nishimori line TO . (This condition of the Nishiniori line will be expained in Sec. 4.) 
It is straightforward to express each contribution to the above partition function Z pictorially in 
terms of the domain walls (see Fig. 6). From the above discussion, it is obvious that if the RBIM is 
in the disordered phase, arbitrary large E' appears and successful error corrections are impossible. 
On the other hand, we can perform successful error corrections if the RBIM is in the ordered phase. 
The 2D RBIM has been studied well and its phase structure is known. When errors and noises 
throughout the measurement of the check operators are totally ignored, the accuracy threshold of 
the toric quantum memory is given by the intersection point of the phase boundary and the Nishimori 
line of the RBIM which is estimated as Pc ~ 0.109 12 . 

Now let us proceed to the second step by taking into account the errors in measuring the check 
operators themseves. In this case, the error chains E live in the two-dimensional space and the 
one-dimensional time because we must repeat the check operations (syndrome measurements) many 
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times. Boundaries of the error chain E reside on sites of the original 3D lattice. In the 3D cubic 
lattice, a link of the original lattice corresponds to a plaquette of the dual lattice. Then the constraint 
(Pn|l can be solved as 

Ui^Y[ai., (23) 

p* 

where the plaquette P* on the dual lattice corresponds to the link £ on the original lattice and 
(T£. is the Ising variable on the link £* of the dual lattice. From Eas. (|17|l and (|23|l . the conditional 
probability proh{E' \E) can be expressed as 

prob(£;) -proM^'lS) =[]exp(J?7p. [] a,.), (24) 
p* (*^p* 

where rjpt = jy^ = ±1 and other notations are self-evident. Then the statistical- mechanical model in 
the present case is the RPGM defined on the 3D lattice. Its third direction represents the real time 
directions along which measurements are done, while each plane lying in the 1 — 2 directions is the 
dual lattice of the origina 2D lattice. 

In the following section, we shall study the phase structure of the RPGM by numerical methods. 
Hereafter we assume that the error probability of syndrome measurement is also given by p for 
simplicity. We shall consider 313 lattice as large as possible. As discussed in Refs.|310], the phase 
structure of an arbitrary large system gives the accuracy threshold of the toric code for the error 
correction. 

4 Phase structure of the RPGM 

In this section, we study the phase struture of the 3D RPGM by MC simulations, and calculate the 
accuracy threshold for the toric quantum memory. The phase boundary is determined by calculating 
the specific heat and expectation values of the Wilson loop. We note that our choice of the above 
model implies that we take into account possible errors which may occur in the measurements of 
the check operators. 

The partition function of the 3D RPGM is given by 

{a^i = ±\\ P P 

where E is the energy and a^i — ±1 is a Z2 classical variable sitting on the link (x, i) of the 3D dual 
lattice (i — 1, 2, 3).^ Hp denotes the product of four cr's along the four links forming P. -qp = ±1 

^As discussed in the previous sections, the present RPGM is defined on the dual lattice. In this section however, 
we shall use simple notations like P instead of P* for the dual plaquette, etc. 
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is a quenched random variable on the plaquette P. The concentration of "wrong-sign" plaquettes 
with i]p = —1 is p as in the previous sections. Thus the weight function P{Kp,ri) for a set of jyp's 
is given by 

P{Kp,v) - l[[{l-p) 5^,s+p6^,^-i] 
p 

^ {2co8hKp)-'^'' -expl^Kp^rjpy (26) 

p 

where Np is the total number of plaquettes, and the parameter Kp is related with p through the 
relations, 

l-poc exp{Kp), p oc exp{~Kp), = e^^^. (27) 

P 

Then the ensemble-averaged free energy F{/3, Kp) is defined as 

F{P,Kp) = -^ PiKp,T^) \ii Zip, r,). (28) 

We presented in Fig.l a schematic phase diagram of the RPGM in the p — T plane. At p = 0, i.e., 
the Z2 gauge system in 3D with a constant gauge coupling, it is known that there are two phases, 
one is the Higgs phase and the other is the confinement phase The critical coupling is estimated 
as (3c = 0.76(Tc = 1.31). As the concentration p of "wrong-sign" plaquettes increases, the critical 
temperature Tc{p) at p decreases and at a certain value of p — po, the critical temperature vanishes, 
Tc — 0. As discussed in the previous section, the system must be in the ordered Higgs phase for the 
toric code to work as a reliable quantum memory. In Fig.l the Nishimori line T = rN(p) is plotted 
by the dashed line. In the high-T region above Tn(p), the critical line is rather horizontal, whereas 
in the low-T region below T^{p), it is almost vertical. 

The Nishimori line is defined in the p — T plane by setting Kp — /3, on which the quantum fluctu- 
ations of gauge field a^i and fluctuations caused by the "random impurity" r]p valance in magnitude, 
i.e., giving the same effect to disorder the system. Let us explain this briefly. For the case p ~ 0, 
there is no randomness, and the quantum(thermal) fluctutions are contoroUed by the Botzmann 
factor exp(/3 (t). On the other hand at T with flnite Kp, the quantum fluctuations are 

suppressed; axi sue frozen and synchronized as Y[p <^ — VP due to the factor expiPrip Y[p c) in Z. 
Thus the effects of randomness is controlled by the factor P{Kp, rj) cx exp(A'p J2p Yi v)- These two 
factors become equal at A'p — (3. 

Let us choose the statistical-mechanical models as the method of error corrections in quantum 
computations. As discussed in Sec. 3, we generate error-correcting chains E' exactly in the same way 
with error chains E, and this procedure corresponds to the Nishimori line in the p — T plane. The 
critical line and the Nishimori line intersect with each other at {pc, T^{pc)). In order that the toric 
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quantum codes to be fault-torelant, the parameters of the toric codes, i.e., {p,T), must be in the 
Higgs phase. The intersection point {pc,T^{pc)) determines the maximum possible value of p for 
fault-tolerant codes, the accuracy threshold, as Pc- 

Let us first focus on the high-T region. In Fig. 7, we plot the specific heat per site, c(/3,p), as a 
function of T for various samples with p = 0.005, 0.01, 0.015, 0.02 and 0.025. For p = 0.005 the peak 
of the specific heat develops as we increases the system size just as in the well studied case of p = 0. 
However, for p = 0.01, the peak does not develop as shown in Fig. 8. For larger p's, the peak of the 
specific heat disappears, although there seems to exist some discontinuities in the derivative dc/df3. 
Below we shall examine these anomalous behaviors of c{j3,p) by calculating the expectation values 
of the Wilson loops, and conclude that there certainly exist phase transitions at these values of p. 
The locations of these phase transitions are pointed with the arrows in Fig. 7. The order of these 
transitions are weak second order or may be higher than 2. In Fig. 9 we plot the critical temperature 
Tc(p) of these phase transitions in the p — T plane together with the solid Nishimori line. Tciji) 
starts from Tc{{)) — 1/0.76 — 1.31 aX p = Q and decreases as p increases, as we expected. 

To confirm the existence of the above phase transitions, we calculated the expectation values of 
the Wilson loop, 

{»?p} {o-xi} C 

= [W'.ample(C^)]cns, (29) 

W^sample(C)= ^ / Z{P.ll), (30) 

{<T„} C 

where C denotes a closed loop in the iD lattice and Jlc ''^a:* product of axiS along C. 

Wsample(C') is the expectation value of the Wilson loop for each sample with fixed {rjp} and [OJcns 
denotes the ensemble average over samples of the quantity O calculated for each sample. W{C) is 
introduced by Wilson to characterize the confinement phase by the so called area law|H], that is, 
W{C) behaves for large C's as 

W{C) oc exp(-aA(C)), (31) 

where A{C) is the smallest area of all the branes whose boundaries are C, and a is a constant called 
the string tension. On the other hand, in the deconfinement phase like the Higgs phase, W{C) 
exhibits the following perimeter law for large C's, 

W{C) cx exp(-7P(C)), (32) 

where P{C) is the length of C. In Fig.lO W{C) are plotted for p ^ 0,0.01 and 0.02 from the 
above. We chose a pair of /3's for each p, slightly above and below Tc{p) that are determined by 
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the calculation of the specific heat (see Fig. 7). In the left column of Fig. 10, W{C) are plotted as 
— \n{W{C))/A{C) vs A{C), so the curves become constant a when the area law (|31|l holds. In 
the right column of Fig. 10, W{C) are plotted as —\n{W{C))/P{C) vs P{C), so the curves become 
constant 7 when the perimeter law H32I) holds. At all the pairs of /3's presented in the figures, we 
observe the changes in the behavior of W{C) from the perimeter law at T < Tc{p) to the area 
law at Tc{p) < T. These changes take place at temperatures very close to Tc, which confirm the 
existence of the phase transitions. The obtained value of Tc{p) by the Wilson loop is insensitive to 
the distribution of the wrong-sign plaquettes, i.e., {?7p} for fixed p. 

Next we study the low-T region. We calculated the specific heat per site c{f3,p) along a line of 
fixed T by changing p, i.e., concentration of the wrong-sign plaquettes. For each p, we prepared 
1000 samples, each of which has a definite configuration of {r/p}, and calculated the specific heat 
for each sample and averaged over the results. To obtain the suitable initial configuration of {cr^i} 
for each sample at low T, we employed the quenched- annealing method. In Fig.ll, the averaged 
values of the specific heat for /3 = 2.5 are plotted from p ~ 0.005 to p = 0.06. We superimpose the 
fiuctuations of each sample around the average. There exists a samll discontinuity in c(/3 = 2.5, p) 
a.t p = 0.027 ^ 0.028. As p increases, the fuctuations also reduce suddenly at p = 0.027 ~ 0.028 
and they become almost constant for p > 0.032 ^ 0.033, so some kind of transition seems to occur 
there. As we show shortly, the same behavior of the specific heat is observed &t (3 — 2.0. In order 
to verify the above expectation, we calculated the expectation values of the Wilson loop, W{C) of 
(I29|l . More precisely, we first calculated expectation value of the Wilson loop W^samplolC*) for each 
sample with fixed {jyp}. Typical results are shown in Fig. 12 for p = 0.013, 0.024, 0.025 and 0.030. It 
is obvious that a change of behavior of W^samplo(C') takes place from the perimeter law (|32|) to the 
area law ()31|l as p increases as we expected. In the region close to p — 0.026, some samples show the 
area law, whereas others exhibit the perimeter law. This result is in contrast to that of the high-T 
region. 

In order to understand the implication of the above result, let us see the expectation value of 
the Wilson loop more precisely. From Ea. (|29|l . it is obvious that the samples with the perimeter 
law dominate over those with the area low as e^^*^*^^ ^ e^'^''-^^ for large C and W{C) exhibits the 
perimeter law if samples with the perimeter and area laws coexist. We then examined Wsample(C') 
for various but a finite number of samples by varying p and found that certain samples at p = 0.025 
and 0.026 exhibit the perimeter law and others show the area law, i.e., the coexistence (we show 
the result of a sample with the area law for p = 0.025 in Fig. 12). For larger values of p, we found 
only the area-law decaying Wsample(C'). (This does not mean that all samples exhibit the area law 
because we only examined a finite number of samples.) In this way we estimate a lower bound of the 
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critical value of p as p ^ 0.026, which is consistent with the calculations of the specific heat given in 
Fig. 11. 

We also studied the case of /3 = 2.0 which is close to the Nishimori line. In Fig. 13, the specific 
heat and the fluctuations over samples are shown. On the Nishimori line, it is proved^ni ^^^X the 
specific heat exhibits no singular behaviors (no divergences). In fact. Fig. 13 shows rather smooth 
behavior in c{j3,p). Let us examine the possibility of a phase transition by studying W{C). In Fig. 14 
we present W{C) at /3 = 2.0 for p = 0.020,0.030,0.031 and 0.040. As p increases, we observe that 
the behavior of W{C) changes from the perimeter to area law. As far as we observed, the maximum 
value of p at which W{C) exhibits the perimeter law is 0.030, and therefore we conclude that the 
phase transition from the Higgs phase to the confinement phase occurs near this value. 

From the calculations of the specific heat and its fluctuations shown in Figs. 11 and 13, we 
estimate the critical values of p as p = 0.032 ^ 0.033 at both (3 = 2.5 and (3 = 2.0. In both cases, the 
fluctuations of the speciflc heat become constant for p > 0.032 ~ 0.033. For the case oi f3 — 2.5, if 
we regard the samll discontinuity in the speciflc heat as the phase transition point, the critical value 
of p is estimated as ^ 0.027. In order to obtain a deflnitive value of Pc^ a thorough investigation on 
the Wilson loop is required. 

5 Conclusion 

In this paper, we studied the RPGM in 3-D numerically and obtained the phase transition line in 
the p — T plane which separates the Higgs and conflnement phases. The critical concentration of the 
"wrong-sign" plaquettes is estimated a.s p — 0.032 ~ 0.033 at low T and also close to the Nishimori 
line. This result determines the accuracy threshold for the quantum toric code as Pc — 0.033. In 
Ref.0], the critical value of p along the T — axis was estimated as 0.029, which corresponds pa in 
Fig.l. This value gives a lower bound of pc, and Fig. 9 shows that our MC simulations are consistent 
with their calculations. 

In Ref. 7, a qualitative phase diagram of the 3D Z2 RPGM was obtained by using the replica 
methods, and the existence of the "gauge-glass phase" is predicted at relatively low T and moderate 
value of p. In the present calculations, however, there appears no sign of that phase. In the higher- 
dimensional gauge models, the gauge-glass phase may exist. 

It was recently suggested that anyonic excitations in Kitaev's model with discrete non-Abelian 
gauge groups may play an important role in constructing quantum gates of QC'sjHlEl- Then gener- 
alization of the present study of the Z2 model to a non-Abelian gauge model is interesting. 
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Figure 1: A schematic phase diagram of the 3D Zi RPGM in the p — T plane, where p is the 
concentration of wrong-sign plaquettes and T = is the "temperature" (the gauge coupHng 
constant). The thick curve is the phase-transition hne Tc{p) which separates the ordered Higgs 
phase and the disordered confinement phase. The dashed Hne Tim(p) is the Nishimori hne which 
corresponds to the present methods of the error corrections in quantum computations. The value Pc 
at the intersection point {pc,T^{pc)) determines the accuracy threshold for the quantum toric code. 
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Figure 2: The original lattice (a) with periodic boundary condition and its dual lattice (b). C^, C\, 
arc two noncontractible closed loops on (a) in the horizontal and vertical directions, while C\ , C\ 
are corresponding loops on the dual site. 
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Figure 3: Path Cx^y* connecting dual sites x* and y* . 
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Figure 4: Typical example of error syndrome produced by noises. The error sites are marked by 
filled circles, where the check operators arc = —1. They are connected by the error links marked 
by thick links, on which extra Z^i's, are applied. 
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Figure 5: An error chain E produced by noises and an applied chain E' produced in the process of 
the error correction. Both E and E' connect the pair of error sites and E + E' form a closed loop, 
i.e., a cycle. 
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Figure 6: Domain wall is produced along the cycle of the error chain E and the error-correction 
chain E' . Corresponding configuration in the partition function l|22|) is achieved by the configuration 
of the dual variables, ax* = —1 within the wall and ax* = +1 outside the wall. 
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Figure 7: Specific heat per site in the high-T region vs In the MC simulation, the size of the 3D 
lattice is Ny — 24'^ sites, and the numbers of sweeps are 20000 for thermalization and 100000 for 
measurement for p=0.005, 0.010 and 0.015, and 30000+100000 for p=0.020 and 0.025. We averaged 
over 10 40 samples with different configurations of rjp. The positions marked by arrows denote 
the locations of phase transitions confirmed by measuring the Wison loops. 
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Figure 8: Specific heat vs /3 for p = 0.005, O.Of . For p = 0.005, tlie peak develops as the system size 
increases, exhibitting a second-order transition. For p = O.Of, the size dependence is very weak; and 
the order of the transition seems to be very weak second order or higher than second-order. 
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Figure 9: Phase diagram in the p — T plane. The crosses denote the phase transition points Tc{p) 
determined by the specfic heat and the Wilson loops. The solid line is the Nishimori line Tn(p). As 
T is lowered, the critical value of p starts to increase, and then decreases. This bending behavior of 
Tc{p) makes the problem of finding the best method of error corrections nontrivial. The statistical- 
mechanical methods of the error correction discussed in Sec. 3 correspond to the Nishimori line, which 
gives a higher accuracy threshold Pc — 0.033 than po — 0.029 determined by Tc{po) = 00. 
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(a) p = 0, Area Law 



(b) p = 0, Perimeter Law 
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(d) p = 0.010, Perimeter Law 
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(e) p = 0.020, Area Law (f) p = 0.020, Perimeter Law 

Figure 10: Wilson loops W{C) in high-T region. The size of the lattice is Ny — 2A^. The results are 
plotted in the pair of axes, A(C)(the smallest number of the plaquettes in all the branes bounded 
by the loop C) and — \n{W (C)) / A{C) in the left column, and P(C)(the number of links contained 
in C) and - ln(W^(C))/P(C) in the right column. If the area(perimeter) law (EIKlISll)) holds, the 
curves in the left(right) column become constants a{"f). 
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Figure 11: Specific heat in the low-T region at (3 = 2.5 vs p. The histogram denotes the fluctuations 
at eachp over 1000 samples around the averge. As p increases, the fluctuation decreases and suddenly 
reduces to an almost constant value at p ~ 0.033 marked by arrows. This value almost coincides 
with the critical value determined by W{C) in Fig. 12. 
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Figure 12: Wilson loop Wsample(C) for certain samples in the low-T region at (3 = 2.5. The size of 
the lattice is Ny = 16^. At p = 0.025 it exhibits the area law, while at p = 0.024, the perimeter law. 
This result gives a lower-bound estimate of the value of p on criticality, which is consistent with the 
speciflc-heat calculations in Fig. 11. (See text.) 
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Figure 13: Specific heat and fluctuations a.t f3 — 2.0. Similar plot as in Fig. 11. As p increases, the 
fuctuation reduces gradually, and reaches at a constant value at p ~ 0.033 as marked with arrows. 
This value almost coincides with the critical value studied by the Wilson loops in Fig. 14. 



(3 = 2.0 



p=0.020 
p=0.030 
p=0.031 
p=0.040 



C 0.02 - 
' 0.015 - 



A:Aroa 

(a) /3 = 2.0, Area Law 



(b) t3 = 2.0, Perimeter Law 



Figure 14: Wilson loop in the low-T region at /3 = 2.0. The size of the lattice is — 12"^. Simiar 
plot as in Fig. 12. We observe that the change in the behavior of W{C) from the area law to the 
perimeter law takes place at p ~ 0.031. This critical value almost coincides with p ~ 0.033 discussed 
in Fig.l3. 
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